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The many-body localization transition (MBLT) between ergodic and many-body localized phase 
in disordered interacting systems is a subject of much recent interest. Statistics of eigenenergies 
is known to be a powerful probe of crossovers between ergodic and integrable systems in simpler 
examples of quantum chaos. We consider the evolution of the spectral statistics across the MBLT, 
starting with mapping to a Brownian motion process that analytically relates the spectral properties 
to the statistics of matrix elements. We demonstrate that the flow from Wigner-Dyson to Poisson 
statistics is a two-stage process. First, fractal enhancement of matrix elements upon approaching 
the MBLT from the metallic side produces an effective power-law interaction between energy levels, 
and leads to a plasma model for level statistics. At the second stage, the gas of eigenvalues has local 
interaction and level statistics belongs to a semi-Poisson universality class. We verify our findings 
numerically on the XXZ spin chain. We provide a microscopic understanding of the level statistics 
across the MBLT and discuss implications for the transition that are strong constraints on possible 
theories. 
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Introduction. Quantum and statistical mechanics rep¬ 
resent two seemingly rather different approaches to the 
description of complex physical systems. Yet these two 
viewpoints agree for a wide class of isolated quantum 
systems, which are said to thermalize [1, 2]. Determin¬ 
ing the circumstances under which an isolated quantum 
many-body system becomes its own thermal bath and 
thermalizes itself, just as Baron Munchausen could pull 
himself out of a mire by his own hair, perhaps using some 
kind of fluctuation, is an open question. 

Phenomena similar to the emergence of thermaliza- 
tion also occur in few-body quantum systems, which 
frequently show the emergence of so-called quantum 
chaos [3]. There, upon changing parameters/number of 
degrees of freedom, the classical system can go from reg¬ 
ular to chaotic behavior. On a quantum level this re¬ 
sults in changes of level statistics, which has proven to 
be a powerful probe of the system properties in the con¬ 
text of quantum chaos. In particular, there exist two 
standard universal limits: Poisson statistics (PS) and 
Wigner-Dyson level statistics (WDS) [4]. For few-body 
systems, PS applies to systems which are classically inte¬ 
grable and do not have any level repulsion. WDS stems 
from random-matrix theory and holds for generic chaotic 
systems, where energy levels repel each other (i.e., the en¬ 
ergy difference between neighboring levels is statistically 
unlikely to be small compared to the mean level spacing). 

Integrable (non-chaotic) behavior is abundant in the 
context of few-body physics. On the other hand, in the 
many-body world the only non-thermalizing phase (in 
the sense of stability to small perturbations) is repre¬ 
sented by many-body localized (MBL) systems [5, 6]. Re¬ 
cent progress established that thermalization fails in the 
MBL phase due to the existence of extensively many con¬ 
served quantities [7-10]. On the other hand, it is known 


that one can tune the system through a phase transition 
into a thermalizing ergodic phase [11-20]. Below we aim 
to understand the evolution of the level statistics across 
the MBL-to-ergodic transition, gaining insights into the 
breakdown of thermalization. 

Crossover between PS and WD statistics has been 
studied extensively in a single-particle physics context: 
for quantum kicked rotor [21], integrability breaking per¬ 
turbations [22, 23], and single-particle Anderson localiza¬ 
tion transition (ALT) [24-26]. In the many-body prob¬ 
lems, PS to WD crossover is also known to occur upon 
breaking of (quantum) integrability [27]. In most of 
the examples, the PS and WDS are the only two sta¬ 
ble points. The only known exception is the ALT, where 
universal statistics different from PS and WDS emerges 
at the mobility edge [24]. 

The spectral statistics in the case of MBL transition 
was demonstrated to evolve from WDS to PS as one local¬ 
izes the system [11, 28-30], however not much is known 
about the intermediate statistics. The common probe 
used to characterize level statistics across MBLT is an av¬ 
erage ratio of the consecutive energy spacings [11-13, 18]. 
However, this is a single parameter and it does not pro¬ 
vide much insight into the intermediate form of the level 
statistic, nor into physical details of its crossover. 

In this paper we study how the spectral statis¬ 
tics changes across the MBL-delocalization transition. 
In order to build a microscopic understanding of the 
level statistics we generalize Dyson’s Brownian motion 
model ]31], previously applied to the ALT [32], to the 
many-body case. From the mapping to Brownian motion, 
we obtain non-trivial relations between fractality [17- 
20], spectral statistics, and properties of matrix elements 
across the MBLT ]20, 33]. While many features can be 
simultaneously explained in this analysis, one surprise is 
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that there appear to be two different regimes of interme¬ 
diate spectral statistics: in one, the effective interaction 
between energy levels in the plasma model has a variable 
power-law, while in the other, the effective interaction is 
short-ranged but over a variable number of levels. 

Within the picture of Brownian motion [31, 32], the 
level statistics is controlled by the effective interaction 
between energy levels, see Fig. 1. In particular, deep 
in the metal phase, the WD statistics emerges from 
the partition function of a one-dimensional Coulomb 
gas, where particles interact with a logarithmic poten¬ 
tial U{s) = logjsj. At a first stage, upon approaching 
the MBL transition, the effective interaction starts to 
decay as a power-law: U{si — Sj) = jsi — when 

|si ~ S 2 I > A^erg- The power-law interaction changes 
tails of the level statistics, so it can be approximately 
described by the plasma model, and is intermediate be¬ 
tween PS and WDS case. At the second stage, when 
exponent 7 becomes bigger than one, the interaction be¬ 
comes effectively short-ranged, and level spacing distri¬ 
bution tends to the semi-Poisson distribution [34]. In this 
regime it is the range of the interaction which changes 
with disorder/system size. As soon as the range of inter¬ 
actions reaches zero, we arrive at Poisson statistics. 

Before discussing implications of the above picture of 
the level statistics, we justify the proposed cartoon us¬ 
ing both analytic and numeric arguments. In particular, 
we argue that the parameter 7 introduced above can be 
extracted from the properties of the many-body matrix 
elements which decay as a power-law with energy separa¬ 
tion between eigenstates, where 7 < I is the same power 
which controls level statistics. The power-law behavior of 
matrix elements can be viewed as a generalization of the 
Chalker-Daniell scaling of wave function overlap [35] to 
the many-body case, and it is consistent with fractality 
of wave functions near MBLT [17-20]. 

Plasma model for level correlations. In the random 
matrix theory, the joint probability density for random 
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FIG. 1. (top) Random walk in a space of Hamiltonians in¬ 
duces a stochastic process on the eigenenergies. The interac¬ 
tion between eigenlevels is set by a potential energy U{si — Sj). 
(bottom) Evolution of the interaction between levels U{s) 
across the MBL transition determines the level statistics. 


matrix ensembles reads 

p-PH 

= ^^H = Y, W{s.) + U{-% - s,),(I) 

i i<j 

where /? = I for orthogonal matrix ensemble which will 
be of primary interest. The confining potential W{s) = 
s ^/2 is parabolic, and interaction is U{si — Sj) = — In js — 
Sjj. As Dyson demonstrated in his pioneering work [31], 
this distribution function may be viewed as a stationary 
distribution of the stochastic random walk in a space of 
matrices (Hamiltonians). 

To derive the joint distribution of eigenenergies from 
a random walk, one can start from the eigenbasis and 
perform a stochastic step in the space of Hamiltonians, 
induced by AH. Then, we get the energy correction in a 
form 

As„ = + V Vmn = {in\AH\n), (2) 

/ 

m^n 

which is the shift of eigenenergies induced by the per¬ 
turbation AH up to second order. For Gaussian ensem¬ 
bles of random matrices, using {VnmVmn) = f At and 
{Vnnymm) = ^mnAr One Can derive Fokker-Planck equa¬ 
tion (see Supplemental Material [36] for more details). 
Its stationary (equilibrium) solution is given by Eq. (I) 
with logarithmic interaction. 

Dyson’s mapping was generalized to the case of dis¬ 
ordered problems [32]. For such problems, it is natural 
to perform a random walk (RW) in a space of Hamil¬ 
tonians by changing realizations of disorder. As we are 
going to concentrate on properties of a spin chain in a 
random magnetic field, which is coupled to the z com¬ 
ponent of a spin S'f, we take AH = hi{T)Sf, with 
{hi{T)hj{T')) = v‘^S{t — T')Sij. Similar to the case of 
random matrices [3, 31, 36], the two correlators which 
determine the level dynamics are: 

{VnnVmm) = 5d„m = {n\Sf\n){m\Sf\m), (3) 

{VnmVmn) = Sc^m = \ {m\Sf \n)\'^, ( 4 ) 

where we assumed that v'^ = 5/L^ where 6 is the many- 
body level spacing, so that represent unfolded energy 
spectrum. The correlator (3) sets the spectrum of a ran¬ 
dom noise, while spectral function Cnm determines the 
interaction between levels in the ensemble. 

Effective interaction between levels. The RW process 
depends crucially on two correlators Eqs. (3)-(4). To 
make analytic progress we use a mean-field like approxi¬ 
mation [32], assuming that dnm and c„m can be replaced 
by their ensemble averages, 

C{uj) = {CnmS{Sn- Sm-Uj)), ( 5 ) 

(and similar expression for dnm) which now depend only 
on the energy difference between eigenstates. Eor the 
single-particle Anderson localization, the Cnm and dnm 
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necessarily coincide with the wave functions overlaps [32], 
Cnm = dnm oc / dx] (t, x) ^ |'0rn (t, x) ^ . The fractality 
of the wave function near the mobility edge results in 
a power-law enhancement of c{uj) oc AjuH [35, 37]. In 
the case of ALT this enhancement arises because the en¬ 
velope of wave functions nearby in energy lives on the 
same multifractal domain [37]. In the many-body case 
similar enhancement can arise from the fractal structure 
of the wave function in the Hilbert space in a vicinity of 
MELT [17-20]. 

Inspired by the approach recently proposed in Ref. [20], 
we apply the fractal scaling to the matrix elements of a 
local operators. In particular, we assume that the in¬ 
verse participation ratio (IPR), Ii = V J2j V~‘^^, 

where ^2 is generalized fractal dimension, and V = 
exp(sL) is the number of states in the Hilbert space. Us¬ 
ing scaling, we translate the IPR into the scaling with the 
distance in the Hilbert space as oc 

where TZ = exp(sdi^k) grows exponentially with (hum¬ 
ming) distance in the Hilbert space, di^k- From here, 
expressing V via frequency, as dy = d/V = w, we get: 

Ek-uj)) oc Finally, omitting 

the diagonal matrix element, we arrive to the scaling: 

c(a;) oc , 7 =I-d 2 - ( 6 ) 

Note, that we did not discuss the microscopic nature 
of a fractal behavior, although Griffiths (rare-region) ef¬ 
fects [18] in vicinity of MBL transition is one possible 
microscopic scenario. Also, relating d 2 to the properties 
of matrix elements, i.e. exponent n in the scaling [20, 33], 
|Unm| oc exp(—(s -|- k)L) is an interesting question. 

The correlation between diagonal matrix elements, the 
function dnm also shows a power-law dependence. How¬ 
ever, there is an enhancement of dnm for n = m, allowing 
to approximate d(a;) as a delta-function, see SM for ad¬ 
ditional discussion [36]. 

Implications for spectral statistics. Using power-law 
form of c(a;) Eq. ( 6 ), and the delta-function form of d(L>j) 
we can map our model onto the plasma model for the 
level statistics [38], provided 7 < 1. The plasma model 
assumes a power-law interaction potential U{s) = A/|s|^ 
in the joint distribution function (1). It predicts the tails 
of the level statistics P[s) cx exp(—for s » I, 
and variance of the number of levels in a box of size 
N becomes vai N cx iV^, which is intermediate between 
WD-like rigidity var cx log and Poisson case [3, 4]. 

For larger values of 7 > 1 the effective interaction in 
the gas of eigenvalues becomes short range, and mapping 
to the plasma model no longer works. Instead, spectral 
properties now are expected to be well-described by a 
family of semi-Poisson distributions [34], which arise from 
a gas of eigenvalues with a finite-range interaction. They 
predict Poisson-like behavior of the tails of P{s) and level 
compressibility P{s) cx s/ 5 g-(/ 3 ?i+i)s^ vai N cx yA^ 
with X < 1, where h is the range of interactions. Such 


level statistics has been dubbed “critical” in the litera¬ 
ture [39-41] and is believed to describe the level statistics 
at the ALT [25, 26]. 

Using the above intuition, we propose the following 
form of the level spacing distribution and spectral rigidity 
to interpolate between WDS and PS, 

P{s;P,^p) = exp (—C' 2 X^“'’'^) , varA^ = 

( 7 ) 

where the parameter 1 > 7p,7var > 0 controls the tails 
of the statistics and level rigidity, and 1 > /3 > 0 de¬ 
termines the level repulsion. The constants Ci ^2 can be 
fixed by requiring that (1) = (s) = 1. When yp = 0, this 
distribution becomes WD. In the opposite limit, yp —>■ 1, 
distribution (7) becomes a semi-Poission with generic /3. 
For the spectral rigidity our interpolating function also 
can describe the (semi-)Poisson limit, however failing to 
capture logarithmic growth of var N in the WD case. 

Numerical results. We use the XXZ spin chain in a 
random field as a specific model with a previously located 
MBL transition [II] to test our picture of level statistics. 
The Hamiltonian is 

i7xxz = ^S,.S, +^u;,5f, = (8) 

(b> * 

where disorder enters via random fields Wi uniformly dis¬ 
tributed in the interval [—IU;IU]. We perform exact di- 
agonalization for chains of size L = 12,... 16(18) with 
periodic boundary conditions to extract properties of ma¬ 
trix elements (spectral statistics). We use central part of 
the many-body spectrum, which corresponds to energy 
density e = {E - Amin)/(F^max - -U^in) = 0.45 ±0.1. The 
MBL transition at this energy density is believed to occur 
near Wc « 3.6 [13]. To unfold levels, we fit the staircase 
function with a 3rd order polynomial. We use both local 
and global level unfolding schemes [42]. 

We start by discussing the numerical results for aver¬ 
aged c{u}), presented in Fig. 2(a). Upon increasing dis- 



FIG. 2. Averaged function c{E) evolves from being almost 
flat at low disorder (W = 0.5) to a power-law decay. Note 
that for the intermediate values of disorder, matrix element 
is enhanced at small energy difference compared to the limit 
of weak disorder. 
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FIG. 3. (a) Evolution of level spacing distributions as system is tuned towards MBL phase. Points represent data, while solid 
lines are best fits with a two-parameter distribution (7). Red and black dashed lines correspond to Poisson and Wigner-Dyson 
distribution, (b) The exponent 7 p, controlling tails of level statistics, flows with L for W < 2.5, but is constant in vicinity of 
MBLT Wc ~ 3.6. (c) In contrast, f5 controlling the level repulsion, remains constant for W <2, and starts to flow closer to the 
MBLT. 


order, we see the crossover of c(uj) from a constant to a 
power-law decay. As one may expect, this crossover hap¬ 
pens at some scale, JVerg, so that c(a; < JVerg) oc const, 
and decays as a power-law beyond oj > A"erg- The ad¬ 
ditional scale JVerg has a meaning similar to the correla¬ 
tion length, over which ergodicity holds. As A^erg —f 0, 
interaction between levels becomes critical even for the 
smallest separations. 

From the power-law form of c(a;), we expect that level 
spacing distribution for the XXZ spin chain to be well 
described by Eq. (7). Fig. 3(a) illustrates that the flow 
of the level statistics is indeed well captured by Eq. (7). 
The P(s) for disorder FF < 2 is not shown, as it looks 
very similar to WD distribution; since P(s) is influenced 
the most by the interaction between close levels, JVgrg 
must become close to zero before we see the flow in the 
level statistics. In contrast to the level statistics, which 
is influenced by a non-critical part of c(a;), the spectral 
rigidity is expected to be less sensitive to the behavior 
of c(w) at small oj. In SM [36] we show that vai N be¬ 
haves as a power-law (7), and becomes linear for W >2. 
Also, we test that different estimates for exponent 7 show 
reasonable agreement as follows from plasma model. 

Einally, we consider the flow of parameters jp and /? 
with increasing system size, presented in Eig. 3(b)-(c). 
While 7 p controlling tails of the level statistics has a 
strong flow at disorder W < 2.5, at larger disorders yp is 
very close to one and changes little with L. This further 
supports the conclusion that for W > 2.5 the effective 
interaction between energy levels becomes short-ranged 
for the largest accessible system sizes. Consistent with 
our expectation, /3 shown in Eig. 3(c) changes weakly 
when statistics is described by plasma model {W < 2), 
and begins to flow once level interactions are local. 

Discussion and open questions. Using analytical and 
numerical arguments we described the spectral proper¬ 
ties across MBL transition using a two-stage flow picture. 
Note that we need at least two parameters, 7 and A^erg, to 
describe level statistics. This is not surprising if we recall 


that even the case of ALT, the existence of multifractal- 
ity means that to describe the universal properties one 
requires more information beyond the small number of 
critical indices needed for a simple thermodynamic phase 
transition [25, 26]. Below we discuss the implications of 
the proposed picture of the spectral statistics flow. 

At the first stage the “correlation length” Nerg shrinks 
to zero, but the exponent responsible for level interac¬ 
tions 7 is smaller than one. Intuitively, the levels beyond 
correlation length become more and more different, cor¬ 
responding to a gradual breakdown of the ETH. Here the 
level statistics can be described by the effective plasma 
model. Although this model was proposed some time 
ago [38], it does not apply in the case of ALT, despite 
the presence of multifractality near single particle mobil¬ 
ity edge. Hence, to the best of our knowledge, the present 
study is the first physical realization of the plasma model. 

The second stage begins at FF > 2.5, when 7 > 1 
so that interactions between levels are local. Although 
we cannot exclude the finite size effects, the numerical 
estimates for the MBL transition at FFc « 3.6 suggest 
that at the MBL transition interactions between levels 
are local. Thus, we conjecture that level statistics near 
and at the MBLT belongs to the same or similar “critical” 
family as the universal statistics at the ALT [39-41]. This 
also naturally explains why the average ratio of the level 
spacing r = min(5„, (5„+i)/max(i5„, 5„+i) at the MBLT, 
widely used in the literature [11-13, 18], is very close to 
the value expected from PS. 

The semi-Poisson level statistics emerges at the same 
value of disorder where the boundary of the Griffiths 
phase was previously identified in the literature [18], 
FF « 2.5 (Refs. [16, 17] report the onset of ergodicity 
breaking at the same location). The existing theories 
of the MBLT [14, 15] predict extensive entanglement 
and subdiffusive transport in the ergodic phase. The 
wide region of critical statistics near transition may be 
a manifestation of finite size effects (system sizes stud¬ 
ied are smaller that diverging correlation length). In- 
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deed the strong overlaps only between adjacent energy 
levels imply logarithmic transport [20], predicted at the 
MBLT [14, 15]. On the other hand, existence of thermo¬ 
dynamically stable Griffiths phase is another intriguing 
possibility. 

In closing, we have found that Dyson’s mapping of 
level statistics to Brownian motion allows one to under¬ 
stand the spectral statistics in the MBL transition at 
least as well as in the ALT for which it was introduced. 
There are basic differences between the two transitions, 
e.g., several quantities which are uniquely defined at the 
ALT allow inequivalent generalizations to the MBLT. 
There are two steps of the spectral statistics flow, one 
with long-range interactions (the plasma model) and one 
with local interactions, and the boundary between the 
two is found numerically to coincide with the onset of 
a Griffiths phase and subdiffusive transport. Since level 
statistics are known to be the simplest universal probe 
of the transition to quantum chaos in simpler problems, 
understanding the origin and universality of the two-step 
plasma model of level statistics is an important challenge 
for theories of the MBLT. 
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Below we present additional details on the derivation 
of the level statistics from the Brownian motion. In par¬ 
ticular, we discuss the approximation for d{u!) used in the 
main text. In the second part we discuss the behavior of 
the spectral rigidity, and compare various estimates for 
exponent 7 . 


Analytic derivation of level statistics from the 
Brownian motion 

Wigner-Dyson statistics 

Let us begin with reproducing the WD statistics from 
the Brownian motion model. We assume the Brownian 
motion in the space of random matrices. 


we find that s„(r) obeys the following Langevin equation: 

-^=ttn[{s}]+Ur), (SIO) 

where drift term is given by Eq. (S 8 ) and white noise is 
specified by 

- t'). (Sll) 

PRM 

Translating the Langevin equation into the Fokker- 
Planck equation governing the evolution of the joint 
probability distribution function P({s„},r), we get: 


dr 





I 92 
Prm ds'^ 




(S12) 


H{t) = Ho+ f dTV{T), 
Jo 


(SI) 


where matrix V (r) satisfies the following properties: 

{Vnm{T)) = -H„m{T), (S2) 

{Vnm(T)Vmn{T')) = 6 {t - t'), (S3) 

{Vnn{T)VmmiT')) = - T')Snm- (S4) 

PRM 

Parameter I3rm specifies the symmetry class, Prm = 1 
for GOE, and /?km = 2 for GUE. 

For convenience, we fix the basis to coincide with the 
(instantaneous) eigenbasis of H{t). We apply the per¬ 
turbation theory to calculate correction to eigenvalues of 
H{t), {s„} induced by the change in the matrix 


/ T+At 

dTV{T). 

Resulting correction to s„ reads: 


As„ — AVnn + 


|A14. 


m^n 




(S5) 


(S 6 ) 


Averaging this equation over V using the fact that = 


Using explicit expression for the drift, Eq. (S 8 ) we see 
that the joint probability distribution 


^({ 4 )=c^ n 


I \0RM I PrM 2 

|Sn~Sm| GXp ^ 


m<.n 

(S13) 

is a stationary solution of the Fokker-Planck equa¬ 
tion (SI2). This distribution function corresponds to 
Wigner-Dyson statistics. It is equivalent to the parti¬ 
tion function of plasma with logarithmic repulsion and a 
parabolic confining potential. 


Plasma model 

After discussing the random matrix example, we move 
on to the many-body case. There the Brownian motion 
is induced by a random walk over different realizations 
of disorder. For a specific case of XXZ spin chain in a 
random magnetic field, we have 

L 

= {hz{T)hj{T)) = v 2 ( 5 ( t - t ')%. 

(S14) 

Repeating the procedure outlined above for the case of 
random matrices, we get the following Langevin equa¬ 
tion: 


-11 — O'"''* 

(As„) = ^„({s})At, 

(S7) 

^=M„[{4T]+en(r), 

(S15) 

Mn[{s}] = S™ + V 
/ 

m^n 

(S 8 ) 

with drift and noise terms reading: 

Pn[{s},T]= c 2^ 2^ _ , 

(S16) 

Using the second moment 


z=i m^n 


2 

PRM 

(S9) 

(C«(T)?m(T')) = y 51 “ P, 

i—1 

(S17) 




FIG. SI. Normalized diagonal correlator {dnm)/{dnn) as a function of |m — n\ for different values of disorder. Although for 
intermediate disorders {dnm) decays as a power-law with |n — m|, the value of is significantly enhanced, compared to the 
power-law tail. 


where the level spacing 6 appeared since we assume that 
s represents unfolded energy spectrum with mean level 
spacing of one. The correlation functions Cnm and dnm 
are defined by the matrix elements of the operator S'f on 
a site i (which couples to the random magnetic field) in 
the instantaneous eigenbasis: 

CnmiT-) = i(n|5f|TO)(m|5'f|n), (S18) 

dlmi^) = ^(n|5f|n)(m|5'f|m). (S19) 

The mean field approximation employed in the main text 
amounts to replacing Cnn^{T) and with their av¬ 

erage over ensemble, which becomes only a function of 
n — m. Fixing = 6jL to cancel the level spacing from 
resulting equations, we get the following Lagnevin dy¬ 
namics: 

^=^,4{s},T]+Ur), (S20) 

with drift and noise terms 

= (S21) 

m^n 

{in{T)im{T')) = d(s„ - Sm)d{T - t '), (S22) 

expressed via rescaled correlators: 

Assuming a delta-function form of d(a;), and a power- 
law ansatz for c(a;), 

c(w) = T-^, 0 < 7 < 1, (S23) 

|a;|' 

the stationary solution of corresponding Fokker-Planck 
equation reproduces the partition function of plasma 
with interaction potential given by: 

= j^- (S24) 


Numerical results for d{u)) 

In order to check the validity of approximating d{u}) 
by a delta-function, we calculate dnm numerically for the 
XXZ spin chain. Figure SI reveals the evolution of av¬ 
erage {dnm) for different system sizes and disorders. At 
disorder W = 0.5, when system is in the ergodic phase, 
we see that dnm when |m — n| > 0 does not depend 
on |to — n|, which is a manifestation of the eigenstate 
thermalization hypothesis(ETH). Nevertheless, values of 
(dnm) for n m are suppressed compared to dnn- Even 
nearby eigenstates can have average spin of different sign, 
this does not contradict to the ETH which rigorously ap¬ 
plies to coarse-grained observables. 

For larger disorder, W > 1.5, the dnm decays approx¬ 
imately as power-law. Nevertheless, the value of dnn is 
still considerably enhanced, hence a proper approxima¬ 
tion for the d{uj) is d{uj) = CiS{uj) + CijoP [scaling of 
the relative weight of delta function in d(a;) in the ther¬ 
modynamic limit is an interesting and open question]. 
Presence of the delta-function contribution in d(a;) is suf¬ 
ficient to make the spectral statistics different compared 
to the case of Anderson transition. In particular, repeat¬ 
ing the mean-field treatment of Ref. [32], we do not get 
constant spectral form-factor K{t,T) in the limit t —?► 0. 
This suggests that level compressibility is vanishing, con¬ 
sistent with plasma model [38]. 


Spectral rigidity and different estimates for 7 across 
the MBL transition 

To probe the spectral rigidity, we study the behavior 
of the variance of number of levels in the box of size N. 
The variance var N as a function of the box size is shown 
in Fig. S2 for L = 16 spins, along with the power-law fits 
of its behavior. The exponent extracted from the fits, 
7 var below is compared to exponents extracted by other 


means. 
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FIG. S2. Variance of the level number in the box of size N 
becomes essentially linear in N for disorder above W > 2.5. 
To minimize the influence of the unfolding procedure, we show 
and use the part where results agree for both local and global 
unfolding. All data is for L = 16 spins. 



FIG. S3. Exponent 7 extracted from c(tu) agrees reasonably 
with the exponent 7 var from the tails of the level statistic. The 
exponent 7 p governing the tails of the P{s) is consistently 
smaller, but agrees with 7 ' = (cilnc(tj)/dlntj)„=o. 

Three exponents, 7 , 7 p, 7 var 5 obtained from matrix el¬ 
ements, level statistics, and level rigidity respectively, 
are expected to coincide in the region of applicability 
of plasma model. In Fig. S3 we plot these exponents 
as a function of disorder for spin chains with L = 16 
spins. We observe the reasonable agreement between 
the values of 7 var and 7 measured from the tails of level 
statistics. On the other hand, yp is consistently smaller, 
and matches much better the 7 ', logarithmic derivative 
of c(w) at w —>■ 0, confirming that P{s) is more strongly 
influenced by a non-universal part in the c(a;). Both 7 
and 7 ' become larger than one at VF « 2.5, indicating 
that for larger disorders the level statistics enters a semi- 
Poisson regime. 










